o 



Exact quantum statistics for electronically nonadiabatic systems using 
continuous path variables 

Nandini Ananth 1 and Thomas F. Miller III 1 

Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, 
USA 

(Dated: 28 July 2011) 

We derive an exact, continuous- variable path integral (PI) representation of the canonical partition function 
for electronically nonadiabatic systems. Utilizing the Stock- Thoss (ST) mapping for an /V-level system, 
matrix elements of the Boltzmann operator are expressed in Cartesian coordinates for both the nuclear and 
electronic degrees of freedom. The PI discretization presented here properly constrains the electronic Cartesian 
coordinates to the physical subspace of the mapping. We numerically demonstrate that the resulting PI-ST 
representation is exact for the calculation of equilibrium properties of systems with coupled electronic and 
nuclear degrees of freedom. We further show that the PI-ST formulation provides a natural means to initialize 
semiclassical trajectories for the calculation of real-time thermal correlation functions, which is numerically 
demonstrated in applications to a series of nonadiabatic model systems. 
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I. INTRODUCTION 

Electronically nonadiabatic processes lie at the heart 
of chemical phenomena, including electron solvation 
dynamicsP^ energy transfer at metal surfaces.^ and ra- 
diationless transitions in the condensed phase!® Eluci- 
dating the mechanisms and rates for these processes re- 
mains a critical challenge from both the theoretical and 
experimental perspectives. 

In an electronically nonadiabatic process, the Born- 
Oppenheimer separation of nuclear and electronic mo- 
tions breaks down, necessitating a description of the cou- 
pled motions of the electrons and nuclei. The exponen- 
tial scaling of exact quantum mechanical methods has 
motivated the development of numerous mixed quantum- 
classical (MQC) methods for nonadiabatic dynamics, in 
which the nuclei are typically treated using classical 
mechanics and electronic degrees of freedom (DoF) are 
treated at the quantum mechanical level. These meth- 
ods, which include the broad classes of mean-fiel d 10 * 11 ! 
and surface hoppingSSHlH approaches, have been succes- 
fully employed in a range of applications. However, it 
has been shown that processes including non-radiative 
electronic relaxation^ and resonance energy transfer^ 
require a consistent description of the coupling between 
the electronic and nuclear DoF, an inherently challenging 
task in the MQC framework. 22 

Semiclassical (SC) methods allow for a dynami- 
cally c onsist ent treatment of the electronic and nuclear 
motionPii^MISl This can be achieved by mapping dis- 
crete electronic states to a continuous variable repre- 
sentation using, for instance, spin coherent stateiP^HH] 
or bos oniza tion techniques based on angular momentum 
theoryP35l] l n particular, Stock and Thoss used the lat- 
ter approach to derive the mapping Hamiltonian, an ex- 
act Cartesian representation of the quantum Hamiltonian 
for an ./V-level system™ The mapping Hamiltonian has 
been successfully used in the SC description of processes 
for which a system initially occupies a pure electronic 
statePSJESl However, use of this approach to calculate 



real-time quantum thermal correlation functions (TCFs) 
has relied on initializing SC trajectories to an approxi- 
mate description of the Boltzmann distributionP3S3 The 
demonstrated sensitivity 25 39 40 of the calculated TCFs to 
the initialization scheme indicates that progress is needed 
to improve the accuracy and generality of this approach. 
In this paper, we derive an exact continuous variable PI 
representation for the Boltzmann distribution of general, 
N-level systems. This is achieved in the mapping frame- 
work using a projection operator to constrain the elec- 
tronic coordinates to the physical subspace of the map- 
ping. We numerically demonstrate that the resulting PI- 
ST representation is exact for the calculation of equilib- 
rium properties of two- and three-state systems with cou- 
pled electronic and nuclear DoF. We further show that 
the PI-ST formulation can be used to initialize SC tra- 
jectories to an exact quantum Boltzmann distribution, 
and using the SC initial value representation (SC-IVR) 
methodpSwe obtain accurate real-time TCFs for a series 
of nonadiabatic model systems. 



II. THEORY 



A. The mapping Hamiltonian 



Consider the general ./V-level Hamiltonian operator 

N 
H = h (R,P) + Y, WWnXVVnl, (1) 



where (R, P) represent the nuclear positions and mo- 
menta, {ip n } is the basis of electronic states, and 
{Kim(R)} is the set of potential energy matrix elements. 
Furthermore, h (R,P) = T(P) + Vo(R), where T(P) is 
the nuclear kinetic energy operator and Vo(R) is a state- 
independent part of the potential energy. 

Following the ST mapping approachj^the ]V-level sys- 
tem is represented by a system of N uncoupled harmonic 



oscillators (HO), such that 



\ipn)(ipm\ -> a+a m 

|^n>->|0 1 2 ---l n ---0 JV ). 



(2) 
(3) 



Here, we have introduced the boson creation and annihi- 
lation operators a+ and a n , which obey the commutation 
rules [a+, a m ] — 5 nm . We have also introduced the singly 
excited oscillator (SEO) states \n) = |0i • • • l n - ■ ■ On), 
which are iV-oscillator eigenstates with a single quantum 
of excitation in the n th mode. The resulting form of the 
Hamiltonian operator is 



H = h (R,P) 



N 

E a+V nrn (R)a m , (4) 

n,m— 1 



or equivalently in the SEO basis, 

N 



H = h (R,P) + E \n)V nm (R)(r 



(5) 



n,m— 1 



Introducing the Cartesian representation of the boson op- 
erators, 

1 i 

x n = -p-{a n + a+) and p n = — = (a+ - a n ), (6) 

we obtain the corresponding Cartesian representation of 
the Hamiltonian operator, 

1 N 

H = h (R, P) + - E {XnX m +PnPm - S nm ) V„(R). 

(7) 
The mapping Hamiltonian in Eq. (17]), also known as 
the Meyer-Miller-Stock-Thoss Hamiltonian, was origi- 
nally derived as a classical model for an electronically 
nonadiabatic system*^ it was later shown to be an exact 
and general representation for the quantum mechanical 
HamiltonianPB 



B. PI discretization 

The canonical partition function is obtained from the 
trace of the Boltzmann operator, 



Z = Tr [e-? H ] 



(8) 



where /3 is the reciprocal temperature and the trace is 
taken over the states that span the electronic and nuclear 
DoF. The resolution of the identity for this space can be 
expressed as 



N 

1= / dR^|R,n}(R,n|, 



(9) 



where R indicates nuclear positions and n indicates the 
SEO state, as before. Repeated insertion of this com- 
pleteness relation yields a path integral discretization of 
the partition function, 



- N P 

I d{R a } E H(R a n a \e-P pH \R c 



Z = (10) 

N P 

la+in a +i) 

{n Q }=la=l 

where P is the number of time-slices and /3p = (3 /P. We 
have introduced the notation J d {R Q } = I Ila=i / dR a ) 

and £{L}=1 = (nLi £n„=i)- 

The standard Trotter approximation^ can be used to 
factorize the matrix elements in this equation, yielding 



P^c 

x exp 

N 



Z= lim /d{R Q }n 



/MP\ 



,-0 P Vo(R a ) 



\\*KP) 



MP 

p 



(Rq — Ra+l) ' (Ra — R«+l) 



X 



e n^i e ^ pv(Ra) i^ + i>, 



(ii) 



{n Q } = l a=l 

where / is the number of nuclear DoF, V(R) = 
£n m=i \ n )Vnm(R)('m\ is the potential energy operator, 
and we set H = 1 throughout this paper. In Eq. ( |TT| ), 
we have assumed that the Hamiltonian in Eq. (fTl) is ex- 
pressed in the diabatic representation, although a sim- 
ilar approach could also be pursued in the adiabatic 
representation.^ Similar expressions for PI discretiza- 
tion in the basis of discrete electronic states have been 
obtainedPEH 



The SEO basis in Eq. (11) can be transformed to the 



Cartesian coordinate basis by first using a projection op- 
erator to select the subset of SEO states from the full set 
of N-oscillator states, 



N 



N 



E 1*1 = 11 



E i-fcW* 

ji=0 



V. 



(12) 



n— 1 i— 1 

where the projectior operator is defined as 

N 

7> = E|n)R 
Then, the full oscillator basis is replaced using 

oo . 

E \3i)Ui\ = / *Ci|xi)(xi|, 

a . — n ** 



(13) 



(14) 



yielding a transformation from the SEO basis to the 
Cartesian coordinate basis for the electronic DoF, 



JV 



Ei* 



dx|x)(x|-p. 



(15) 



As in Klauder's work with spin coherent states}**" the pro- states, 



jection operator in Eq. (15) constrains electronic coordi- 



nates to a specific manifold in phase space. However, PI 
formulations using spin coherent states have only prov en 
numerically tractable in the semiclassical limit j22I4aEl3 
whereas we shall derive a PI formulation that can be 
used for exact numerical simulations. 



Substituting Eq. ( 15 ) into Eq. (Ill, we obtain an exact wn pre 



PI discretization of the canonical partition function in 
continuous variables, 



z-fcoj^iun^y 



X e ~Trtr( R «~ R a+i) -( R o 



-Rc+i) 



o -/3pV (Rc,) 



d{x a } 



n 



I I {Xn 



D -0 P V(R a 



'P|x aH 



(16) 



leaving only the task of evaluating the matrix elements 
in the last term to obtain a computationally useful ex- 
pression. 

We note that a short-time approximation to the elec- 
tronic matrix elements could be performed directly in the 
Cartesian representation at this stage. However, we find 
that a more numerically stable result is obtained by mak- 
ing the approximation in the SEO representation, such 
that 



N 



( x | e -/^V(R).p |x A = 



(x|n)7W nm (R)(m|x'),(17) 



where .M nTO (R) = (n|e~^ pV ( R )|m). Recognizing that 
the coordinate space SEO wavefunction is the product 
of (N — 1) ground state HO wavefunctions and one first 
excited state HO wavefunction, we have 



V2 



(18) 



where [.] n denotes the n th component of the enclosed 
vector. The Boltzmann matrix element in the SEO 
representation can then be obtained following textbook 
procedures p2 such that to order #(/3|>), 



MnmCR) 



-/9 P V„„(R) 



, n = m, 



-P P V nm (R) e-^-W ,n^m. 



(19) 



In the zero coupling limit, the matrix elements .A/f nm (R) 
assume a diagonal form so that the different compo- 
nents of the electronic position vectors do not mix. 
The off-diagonal matrix elements are related to the 
penalty of ring-polymer kink formation, in which neigh- 
boring PI tim e-slic es reside on different diabatic elec- 
tronic surfaces] 52 * 53 ^ 



Finally, substituting Eqs. (18) and (19) into Eq. (16), 



we arrive at the exact, continuous- variable PI-ST repre- 
sentation of the canonical partition function for a nonadi- 
abatic system with / nuclear DoF coupled to N electronic 



lira 



IMP 
P'~ao \/3n N + 1 
p 

I J. a a -^ Q ' 



d{TL a } / d{x Q } 



(20) 



A a =e _ ^r (R "* _:R <*+ l) -( R «- R «+i) p -/3pV (R a ) 



, , i -i n..-a. a+ i) e -ppv (n. a ) (21) 



F a =*a M(R a ) x Q+1 , and 



(22) 
(23) 



Eqs. (20) - (23) will be used to calculate numerically 



exact equilibrium properties for nonadiabatic systems. 

We note that our PI-ST formulation is different from 
the result derived in Ref. 54, since we include a projection 
operator to constrain the system to the physical subspace 
of the mapping. For cases where the system is prepared 
in a pure SEO state, this constraint is implicitly obeyed, 
and the two PI representations are equivalent. However, 
treatment of Boltzmann distributed systems require that 
the projection operator be explicitly included, as is done 
in the present study. 



III. EQUILIBRIUM SIMULATIONS 
A. Implementation details 

The equilibrium properties considered in the current 
paper include the nuclear probability distribution, the 
state-specific nuclear probability distribution, and the 
average total energy. All equilibrium simulations were 
performed using standard path integral Monte Carlo 
(PIMC) importance sampling techniques, although the 
use of path integral molecular dynamics (PIMD) meth- 
ods is also straightforward with the formulation devel- 
oped here. 

The nuclear probability distribution is defined as 



P(R) 



Tr[(S(R - K)e-' 3H } 
TY[e-f> H ] 



(24) 



Using the PI-ST representation for the Boltzmann oper- 
ator, we obtain 



P(R)= 



/ d{R a } J d{x a }5(n - Rp) H A a Q a T a 

q = l 

p ■ 

J d {R a } / {dx Q } J^[ Aa.Ga.Tct 



(25) 



where A a , T a and Q a are defined in Eqs. (21) - (23 1 



Importance sampling can then be performed using 

p 

vf({x q },{r q }) = Y[A a g a \? a \, 



where the absolute value in the last term ensures a non- 
negative sampling function. The expression for the nu- 
clear probability distribution is thus 



P(R) = 



<J(R - Rp)sgn(J-) 



sgn(J r ) 



w 



where 



sgn(J") = Y[ F<x/\ ?<* 



(26) 



(27) 



and the angle brackets in Eq. ( 26 ) indicate the ensemble 



average with respect to the distribution FF({x a }, {R Q }). 
The state- specific nuclear probability distribution is ob- 
tained from the projection of the nuclear probability dis- 
tribution onto a given electronic state, 



P(n,R) = 



Tr[(5(R-R)|n)(n|e-^] 

Tv[e-P H ) 
( 6(R - R P )f n sgn(^) ) w 



sgn(J r ) 



w 



where 



T = 

•J n. 



[x^(Rp)xi], 
x^7W(Rp) Xl 



(28) 



(29) 



The elements of the matrix Ai (R) are defined in Eq. ( 19 1 



The average total energy of the system is obtained us- 
ing a primitive energy estimator, 



(E) 



1 dZ 

~Z~d(3 



%) sgn(-F) 



w 



sgn(J") 



w 



where 



p T -9-M(R„) 



~ ^ xj7W(R Q )x Q+ i 

a — 1 u> \ / \ 



(30) 



(31) 



and 



dA 
<9/? 



E 



— ^-(R Q — Ra+i) '{Ret — Rq+i) 



2f3 2 



V Q (R a ) 



(32) 



B. Equilibrium simulation results 



The first model that we consider (model I) includes a 
two-state system coupled to a single vibrational DoF; it 
is a standard benchmark for the tr eatm ent of equilibrium 
statistics in nonadiabatic systems E2QS 




FIG. 1. Diabatic potential energy curves for model I, with 
state 1 in red (solid line), state 2 in blue (dashed line), and 
the coupling in green (dash-dotted line) 



The matrix elements of the diabatic potential operator 
are 

Vu = -ki(R - Ri) 2 + ei and 

V l3 = ce- a( - R ~ R ^\ (33) 

where the potential energy parameters are specified in 
Table [I] The simulation is performed with a nuclear mass 
of 3600 a.u. and temperature T = 8 K. The potential 
energy curves for Model I are plotted in Fig. [I] 

TABLE I. Parameters for Model I 



Parameter 


Value (in a.u.) 


kx 


4 x 10 -5 


k 2 


3.2 x 1CT 5 


Ri 


-1.75 


R 2 


1.75 


ei 





£2 


2.28 x 10" 5 


c 


5 x 1CT 5 


a 


0.4 


Rl2 






The nuclear probability distribution [Eq. (26)] ob- 



tained from a 32 bead PI-ST simulation is shown in 
Fig. ^ja) and is graphically indistinguishable from a nu- 
merically exact discrete variable representation (DVR) 55 
grid calculation. The tight convergence in this plot was 
achieved using 5 x 10 9 Monte Carlo (MC) steps, and a 
similar number of steps was found to be necessary for the 
corresponding PIMC calculation in the discrete diabatic- 
state representation. In Fig. [2Kb), we show that the 
state-specific nuclear probability distribution from this 
simulation also reproduces the exact results. We further 
calculate the average total energy for Model I and show, 
in Table [ITJ that the PI-ST result approaches the exact 
value in the limit of a large number of beads. 




FIG. 2. (a) The nuclear probability distribution for model 
I, obtained using a 32-bead PI-ST simulation (red, solid line) 
and an exact grid calculation (black squares), (b) The state- 
specific nuclear probability distribution obtained from the PI- 
ST simulation, with state 1 in red (solid line, left peak) and 
state 2 in blue (solid line, right peak); black squares corre- 
spond to the exact results. 




-2 2 

R(a.u.) 



FIG. 3. Diabatic potential energy curves for model II, with 
the G state in red (solid line), the LE state in blue (dashed 
line) and the CT state in green (dash-dotted line); the con- 
stant coupling elements are not shown. 



TABLE III. Parameters for Model II 



Parameter 


Value (in eV) 


OJa 


0.05 


ko 





kcr 


0.1 


kLE 





EG 





ECT 


0.25 


£LE 


0.25 


cg,ct 


0.02 


CCT,LE 


0.03 


CG,LE 






TABLE II. Average Energy for Model I at T = 8 K 



No: of Beads 


Energy (10~ 5 a.u.) 


8 


5.09 


16 


5.12 


32 


5.14 


Exact 


5.145 



a Statistical error for all cases is less than 10 7 a.u. 

Model II is a three-state system coupled to a single 
vibrational DoF. It is based on a model used to simulate 
ultrafast photoinduced electron-transfer P2 The model in- 
cludes a ground (G) electronic state, a locally excited 
(LE) state that is accessible via photoexcitation, and a 
charge transfer (CT) state that facilitates radiationless 
decay to the ground state. The CT state acts as a bridge 
state between the ground and LE states, and it is coupled 
to both these states via a constant potential; there is no 
direct coupling between the ground and LE states. 



The matrix elements of the diabatic potential operator 
are 



1 „ 
Vi = -u s R + kiR + Ei and 

Vij = Cij, 



(34) 



where i,j £ {G,CT,LE}, and the nuclear mass is 
544.23 a.u. The potential energy parameters for Model 
II are provided in Table |III| and the diabatic three-state 
potential is shown in Fig.[3l The simulation is performed 
at T = 1500 K, chosen such that all three states are ther- 
mally accessible. 

The converged nuclear probability distribution for this 
model is obtained from a 4-bead calculation and is graph- 
ically indistinguishable from the exact results from a 
DVR grid calculation, as seen in Fig. |4|a). These results 
were obtained using 10 9 MC steps. The state-specific 
nuclear probability distributions shown in Fig. Elb) also 
reproduce the exact results. 



IV. THERMAL CORRELATION FUNCTIONS 
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FIG. 4. (a) The nuclear probability distribution for model 
II, obtained using a 4-bead PI-ST simulation (red, solid line) 
and an exact grid calculation (black squares), (b) The state- 
specific nuclear probability distribution obtained from the PI- 
ST simulation with the G state, the CT state and the LE 
state in red, green and blue (solid lines, in order of decreasing 
population) respectively. 



Further, in Table [TV] the results of the average energy 
calculation are reported, and the exact results are recov- 
ered with increasing bead numbers. 



TABLE IV. Average Energy for Model II at T= 1500 K 



No: of Beads 


Energy (10~ 3 a.u.) 


1 

2 

4 

Exact 


5.54 
6.63 
6.69 
6.688 



Statistical error in all cases is less than 10 a.u. 



The equilibrium properties calculated for these model 
systems demonstrate that the PI-ST representation pro- 
vides a general and exact statistical description of elec- 
tronically nonadiabatic systems. 



The PI-ST representation provides a natural means to 
initialize SC trajectories to an exact quantum Boltzmann 
distribution for the calculation of real-time TCFs. In this 
paper, we demonstrate this using the SC-IVR method, 
which has already be en succ essfully implemented in 
the mapping framework! 32 ! 35 * 6 ^ However, any trajectory- 
based model for real-time dynamics could be combined 
with our exact PI-ST formulation. 

A general real-time TCF is expressed as 

Cab® = ^Tr [e^ H Ae iHt Be- iHt ] , (35) 

where A and B are generic operators. Substituting the 
PI-ST representation for the Boltzmann operator from 
Eq. pfl, the TCF can be written 



p—i 
Cab® = \ J d{R a } f d{x Q } JJ A a Q a T a 



a=l 
l ' pH AHt-D^-iHt, l3p " 



x (x P ,Rp\e--^- Ae ltlt Be~ %nt e — s-p|xi,Ri), (36) 
where A a , G a , and T a are defined in Eqs. (2Tp3|. 



A. Herman-Kluk (HK) IVR 

The HK-IVR propagatorJSMHH 



,-im 



= (2tt) 



~(N+f) 



dz |z t )Cf K (z )e 



JSt(zo), 



zo|, 



(37) 

is a coherent state approximation to the full coordinate 
state SC-IVR propagator. Here, |z ) = |x ,Po)|Ro, Po) 
represents the initial electronic and nuclear coherent 
states of width 7 and T, respectively, and |z t ) is ob- 
tained from classically time-evolving the initial positions 
and momenta for time t. In addition, St is the classical 
action, and the HK prefactor is given by 



Cf K (z ) = Det 



1 T dz t _ ± 

2 S dz^ 8 



(38) 



whereg=(( 7 ,r) 1 /2 ;l ( 7 ,r)- 1 /2). 



The forward and backward propagators in Eq. ( 36 ) can 



be replaced by HK-IVR propagators to obtain an expres- 
sion with a double phase-space integral over initial con- 
ditions, 



c^(t) 



(27T) 



-2(JV+/) 



f d{R a } I d{ Xa } Y[ A a g a F< 



x / dz J dz' e^ s -^ +s ^C- t (z' )C t (z )(z' t \B\z t ) 

x (xp,R P | e -^ H / 2 ^| z ; ) )(zo| e -^ ff / 2 7'|x 1 ,R 1 ). (39) 

MC integration of the resulting oscillatory integrand is 
known to be challenging;— and despite several advances 



in the evaluation of such integrandsp^nSSl ^he HK-IVR 
approach is limited to systems with few DoF. Nonethe- 
less, we include the HK-IVR implementation to illustrate 
the generality of our exact PI initialization approach and 
to provide a reference semiclassical result to compare 
against the linearized IVR implementation. 



the exact PI-ST representation, we expect that this ap- 
proach would introduce only small deviations from the 
exact quantum statistical description of a nonadiabatic 
system. 



C. Implementation of PI-ST initialization 



B. Linearized IVR (LSC-IVR) 

The LSC-IVR approximation to the coordinate state 
SC-IVR expression for correlation functions is obtained 
from a first order expansion of the difference in the ac- 
tions of the forward and backward trajectories! 64 * 65 * The 
resulting expression corresponds to the classical Wigner 
model and can be written 



c^(ty- 



(27T) 



<N+f) 



(40) 
dp Q J dx Q A^(x ,p )B w (xL t ,Pt), 



where A^ 



(e l3H A), and the Wigner transformed 
operators are obtained by evaluating 



O w (x,p) 



dAx(x 



Ax, 



|0|x- 



Ax, 



3 ip J Ax 



(41) 



The LSC-IVR approximat ion lar gely fails to capture 
quantum coherence effects,^ 41 * 651 but it successfully de- 
scribes other quantum effects such as zero point energy 
and tunneling, m aking it suitable for many condensed 
phase applications! 4116 ^ 1 ^ 

The PI-ST representation of the Boltzmann operator 
can be substituted in the expression for the TCF in 
Eq. J4lj) to obtain 



(o 7r r {N+f) r r p ~ 1 

Cfjf(t) = y -^ / d{R a } I rf{x Q } [] A a G a F a 

a— 1 

x f dz A^(z )B w (z t ), (42) 
where 



I^(z ) = / dAx / dARe ,p °' Ax+,p »' AR 
AR, 



Ax „ 
x (x ^' R o 



|e 2^| XljRl ) 



~ Pp u . Ax „ AR, 

x(xp,R P |e — ^^1x0 + — ,Ro + — ) 



(43) 



We recognize that using the exact PI-ST representa- 
tion of the Boltzmann operator introduces an oscillatory 
term in the LSC-IVR formulation via A^(z ). For fu- 
ture applications to large systems, this oscillatory term 
can be eliminated using techniques such as the Ther- 
mal Gaussian Approximation, for the Boltzmann ma- 
trix elements in Eq. (43). Since the remaining (P — 1) 



Boltzmann terms in Eq. ( 42 ) will still be treated using 



Equations ( 39 1 and ( 42 ) can both be expressed in the 
form 



ci B {t) = 



1 



d{R a } /d{x a }^({x Q },{R Q }) 



x /({x Q },{R a })$«(x 1 ,x P ,R 1 ,Rp,i) 
= ( $«(x 1 ,x P ,R 1 ,Rp,t)/({x Q },{R Q }) ) w 

/(WWNDw. ( 44 ) 
where M^({x a }, {R Q }) is a sampling distribution and 
/({x Q },{R Q }) and /z({x Q }, {R Q }) are weighting fac- 
tors, all of which emerge from the PI-ST treatment of the 
Boltzmann operator. The term $^(x 1; R^xp, Rp, t) in 
Eq. ( 44 ) contains the real-time information obtained from 
the SC trajectories, and the superscript £ € {HK,LSC} 
indicates which SC approximation is employed. 

The calculation of the TCF in Eq. ( 44 ) is performed 



by first generating an ensemble of configurations from 
the probability distribution VF({x Q }, {R Q }). Then, as in 
standard SC-IVR calculations p^ MC importance sam- 
pling is used to evaluate $*(xi,Ri,Xp,Rp,i). For the 
HK-IVR implementation, 

$ HK (x 1 ,xp,R 1 ,R P ,t) = 
7V HK (i)<0 HK (z o ,z i ,z^,z; ; x 1 ,xp,R 1 ,Rp)) I 



(45) 



'n HK , 



and for the LSC-IVR implementation, 



$ 



LSC 



(xi,xp,Ri,R P ,<) = 
Ar LSC (t)(^ LSC ( Z0)Zt ;xi,x F ,Ri,Rp)) 



(46) 



IJLSC 



where n^ is a probability distribution function used to 
generate an ensemble of initial coordinates and momenta 
for the SC trajectories, <jfi is the corresponding time- 
dependent estimator, and N^(i) is the associated time- 
dependent normalization term. 

In this paper we calculate the electronic state popula- 
tion TCF, C|„(t), where A = B = \n){n\ in Eq. ((35 1. 



For this special case, the detailed form of the functions 
described above are provided in the appendix. 



D. Dynamics simulation results 

We calculate the real-time electronic state population 
TCF using the PI-ST representation of the Boltzmann 
distribution to initialize SC trajectories for dynamics in 
the LSC-IVR and HK-IVR frameworks. The first set of 
results presented are for model III, a simple two-state 
system described by the Hamiltonian, 



H 



Aov 



(47) 



where a z and a x are the Pauli matrices. The potential 
parameters are (a, A) = (0.5, 1) in a.u. The mapping 
Hamiltonian for a two-state system assumes a quadratic 
form for which both the HK-IVR and LSC-IVR formula- 
tions are exact. The simulation is performed at a recip- 
rocal temperature of (3 — 1 a.u., using coherent states 
of width 7 = 1 a.u. SC trajectories are integrated 
using an Adams-Bashforth-Moulton predictor-corrector 
integrator.^ Exact results are obtained from a DVR grid 
calculation. 

Fig. [5] illustrates that the Cn(i) TCF calculated using 
the LSC-IVR implementation reproduces the exact re- 
sults, as expected. Simulations performed using the HK- 
IVR implementation yielded graphically indistinguish- 
able results. These calculations were performed using an 
8-bead simulation with 10 8 MC steps for the sampling 
of the probability distribution lV({x, R}). For each of 
these configurations, ensembles of 120 trajectories and 
10 trajectories were generated to obtain the HK-IVR and 
LSC-IVR TCFs, respectively. 



u 




In Fig. [6j the TCFs from the HK-IVR and LSC-IVR 
simulations are compared with the results from an ex- 
act DVR grid calculation. The HK-IVR implementation 
reproduces the exact results with remarkable accuracy, 
although the calculation required an ensemble of 80, 000 
trajectories per equilibrium configuration to converge re- 
sults up to a time of 2.5 a.u.; results for longer times 
would require even larger numbers of trajectories. The 
LSC-IVR simulation, however, required only 8000 tra- 
jectories per equilibrium configuration to achieve conver- 
gence up to 5 a.u. in time. As expected, the LSC-IVR 
approximations dampens the oscillations seen in the ex- 
act calculation. 
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FIG. 6. The real-time electronic state population TCF for 
model IV, obtained from the HK-IVR method with PI-ST 
initialization (red, dotted line), the LSC-IVR method with 
PI-ST initialization (blue, solid line), and an exact grid cal- 
culation (black, dashed line). 



FIG. 5. The real-time electronic state population TCF for 
model III, obtained from the LSC-IVR method with PI-ST 
initialization (red, solid line) and an exact grid calculation 
(black squares). Graphically indistinguishable results were 
also obtained using the HK-IVR method with PI-ST initial- 
ization. 



Model IV is a two-state system coupled to a single 
nuclear DoF of mass M = 1 a.u. The Hamiltonian for 
this model is 



H 



2M 



2 R 



aRa z + Aa x , 



(48) 



with parameters (in a.u.) a = 1, A = 1, and 
f3 = 1. We use coherent states of width 7 = 1 a.u 
and r = 1 a.u. for the electronic and nuclear DoF, re- 
spectively and 4-bead simulations were performed with 
10 8 MC steps for sampling the probability distribution 
W({x,R}). 



V. CONCLUSIONS 

We have derived an exact PI-ST representation for the 
Boltzmann statistics of N-level systems using continuous 
path variables for both the electronic and nuclear DoF. 
This result is demonstrated to be numerically exact for 
equilibrium simulations of two- and three-state systems. 
Additionally, the PI-ST representation is used to initial- 
ize trajectories in the SC-IVR framework allowing for 
the calculation of real-time TCFs with encouraging ac- 
curacy. Natural future applications of this methodology 
include charge transfer reactions in the condensed phase 
and metal-surface energy transfer processes, for which 
excited electronic states are thermally accessible. 
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Appendix A: Electronic state population TCF 

The detailed form for all functions used to calculate 
C|„(4) are provided here. The functions ^({x,,}, {R Q }), 
/({x Q },{R«}), and / z ({x Q }, {R Q }) arise from the PF 
ST representation of the Boltzmann operator and are 
identical for both the LSC-IVR and HK-IVR implemen- 
tations. Specifically, 

W({x a },{R a }) = e -^o(Rp)/2 g p 1^1 

p-1 
X J] e -/»pV (R„) Aa g a pr^ (A1) 

a=l 



In the LSC-IVR framework, the probability distribu- 
tion is 

tLSC 



IT^(zo;Ri,Rp) = 



_ IMP I 

e f> 



:(Ro-i(Rp+Ri)) T -(R -i(Rp+Ri)) 

Xe _ 2MTP P -Po-xJ-Xo-p^-po 

and the corresponding estimator is 



(A6) 



/LSC 



(z ,z t ;xi,x P ,Ri,Rp) = [x t ] + [Pt] n - 



1 



[x + i Po ] n (x - z Po ) T X'(Ri)x 1 - - [M'(Ri)xi] f 



[x T P M'(R P )] i 



iP^(Rp-Ri)-x[-x t - P ;-pt 



(A7) 



A final step in the TCF calculation is to obtain N^ (t) , 
which is required to both normalize the probability distri- 
bution function 11^ and to correct for the non-unitarity of 
the SC propagatorPS Enforcing that the total electronic 
state population is conserved at all times yields 



where A a , Q a , and T a are defined in Eqs. |21p3 l; 

e -/3py (Ri)/2 sgn (^) 



NHt) = 



cL(o) 



2^m=l Cram(i) 



(A8) 



/({x Q },{R Q }) 



x£7W(Rp) Xl 



(A2) 



where sgn(J r ) is defined in Eq.(27l and the elements of 
M defined in Eq. (IT9|; and 



where C| n (0) = Tr \e~^ H V n \ 1% can be calculated from 
an exact PI-S T e quilibrium simulation. The terms 
C^ m (t) in Eq. (A8) are un- normalized TCFs, such that 



-(Rp-Ri) -(Rp-Ri) 



fz ({x Q }, {R«}) = sgn(.T)e 

(A3) 
The terms required to evaluate $£(xi,Ri,xp, Rp,i) 



in Eqs. (46) and (47) are derived by substituting A 
B = |n)(n 



into Eqs. (39 1 and (42). In the HK-IVR 



framework, this yields, the probability distribution func- 
tion 

n HK (z ,z^;R 1 ,Rp)= e -/JrArp( p "- p o+ p o T -Po) 
e ~ /3r+2Afp ((Ro-Rp) T -(Rd-Rp)+(Ro-Ri) T -(Ro-Ri)j 

e ~ 2(7+1) ( x • x Q+ x O T - x o)- 2(-f + l) (P(T-PO+Po T -Po) (^4) 

and the corresponding estimator 

^ HK (z ,Zo,z t ,z' ( ;xi,xp,Ri,Rp) = 
[loA - ip'tln [7o«t + iPt] n 
x [x^M'(Rp)( 7 x^ + zp^)] ri (( 7 ox -*Po) T X'(Ri)xi) 

x C_ i (z^)C , t (z )e 4[s -' (z o) +s *( z «)] 

x e ; ^FT Po ' Xo_ wtPo ' x a _ 4( R 't _R t) -( R t-Rt) 

x e -4L(p;-p t ) r .(p;-p t )+i(p t +p;) r -(R;-R t ) 

x e ^r; A 2Mp (Po T -(Rj-Ro)-p^(Ri-Ro))+^TT(p; T - x ;-pr- x 

X e -2TTFT)( x ^ x « +x t T - x t)~2TTfTT (p ^' Pt+p f T ' p t\ (A5) 

where the elements of the matrix A4'(R) are identical to 
those of the matrix in Eq. ( 19 ) with fi —> (3/2. 



Ci m (t) = N£(i)Ci m (t), where A = \n)(n\ and B = 
\m)(m\ in Eq. ( |35| . These un- normalized terms are ob- 
a pVi ( R toi|ied following the steps described above, except that 
the first line on the right-hand side of Eqs. (A5) and (A7) 
is modified to include the m th (rather than the n th ) com- 
ponent of vectors x t , pt,xj, andpj. The additional com- 
putational cost associated with this normalization term 
is negligible. 
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